Analysis of Tabula Muris with celda

Install appropriate R packages

install.packages(c("renv", "devtools", "BiocManager"))
renv::activate()
renv::restore(lockfile="renv.lock_no_celda")
devtools::install_github("campbio/celda@v1.7.3")
BiocManager::install(c("ExperimentHub", "TabulaMurisData"))

library(SingleCellExperiment)
library(ExperimentHub)
library(TabulaMurisData)
library(celda)

library(Matrix)
eh <- ExperimentHub()

snapshotDate(): 2020-10-27

data <- eh[["EH1617"]]

see ?TabulaMurisData and browseVignettes('TabulaMurisData') for documentation
loading from cache

sce <- data[, data$tissue == "Lung" & is.na(data$subtissue)]
counts(sce) <- as.matrix(counts(sce))

sce <- selectFeatures(sce, minCell = 3, minCount = 3)
moduleSplit <- recursiveSplitModule(sce, sampleLabel = sce$mouse_id, initialL = 10, maxL = 150)

==================================================
Starting recursive module splitting.
==================================================
Sat Jan 29 14:31:36 2022 .. Collapsing to 100 temporary cell populations
Sat Jan 29 14:33:16 2022 .. Initializing with 10 modules
Sat Jan 29 14:33:58 2022 .. Created module 11 | logLik: -81269791.8121076
Sat Jan 29 14:34:08 2022 .. Created module 12 | logLik: -81171328.5303145
Sat Jan 29 14:34:18 2022 .. Created module 13 | logLik: -81104795.1267257
Sat Jan 29 14:34:28 2022 .. Created module 14 | logLik: -81032372.577771
Sat Jan 29 14:34:36 2022 .. Created module 15 | logLik: -81009827.6001511
Sat Jan 29 14:34:43 2022 .. Created module 16 | logLik: -80969733.8288011
Sat Jan 29 14:34:50 2022 .. Created module 17 | logLik: -80944295.1708537
Sat Jan 29 14:34:59 2022 .. Created module 18 | logLik: -80900411.9338321
Sat Jan 29 14:35:09 2022 .. Created module 19 | logLik: -80867745.7719349
Sat Jan 29 14:35:21 2022 .. Created module 20 | logLik: -80832754.0375926
Sat Jan 29 14:35:32 2022 .. Created module 21 | logLik: -80810271.9547412
Sat Jan 29 14:35:44 2022 .. Created module 22 | logLik: -80783992.0275586
Sat Jan 29 14:35:57 2022 .. Created module 23 | logLik: -80763129.1718582
Sat Jan 29 14:36:07 2022 .. Created module 24 | logLik: -80747550.3399068
Sat Jan 29 14:36:15 2022 .. Created module 25 | logLik: -80735226.4809027
Sat Jan 29 14:36:22 2022 .. Created module 26 | logLik: -80723247.5059884
Sat Jan 29 14:36:28 2022 .. Created module 27 | logLik: -80705915.4627029
Sat Jan 29 14:36:39 2022 .. Created module 28 | logLik: -80686644.7080038
Sat Jan 29 14:36:46 2022 .. Created module 29 | logLik: -80677210.3078177
Sat Jan 29 14:36:52 2022 .. Created module 30 | logLik: -80669333.4473654
Sat Jan 29 14:36:59 2022 .. Created module 31 | logLik: -80662182.7162198
Sat Jan 29 14:37:07 2022 .. Created module 32 | logLik: -80651418.2908446
Sat Jan 29 14:37:17 2022 .. Created module 33 | logLik: -80636139.4886655
Sat Jan 29 14:37:27 2022 .. Created module 34 | logLik: -80623405.1537938
Sat Jan 29 14:37:33 2022 .. Created module 35 | logLik: -80618787.0742386
Sat Jan 29 14:37:43 2022 .. Created module 36 | logLik: -80615732.2013244
Sat Jan 29 14:37:55 2022 .. Created module 37 | logLik: -80608020.8869632
Sat Jan 29 14:38:02 2022 .. Created module 38 | logLik: -80607613.1043613
Sat Jan 29 14:38:11 2022 .. Created module 39 | logLik: -80601990.450228
Sat Jan 29 14:38:21 2022 .. Created module 40 | logLik: -80596049.6822151
Sat Jan 29 14:38:29 2022 .. Created module 41 | logLik: -80590385.508503
Sat Jan 29 14:38:36 2022 .. Created module 42 | logLik: -80585143.9416136
Sat Jan 29 14:38:43 2022 .. Created module 43 | logLik: -80581413.2768909
Sat Jan 29 14:38:52 2022 .. Created module 44 | logLik: -80577029.2118674
Sat Jan 29 14:39:01 2022 .. Created module 45 | logLik: -80573818.1305742
Sat Jan 29 14:39:08 2022 .. Created module 46 | logLik: -80570787.615158
Sat Jan 29 14:39:15 2022 .. Created module 47 | logLik: -80565441.1056954
Sat Jan 29 14:39:23 2022 .. Created module 48 | logLik: -80559677.4315346
Sat Jan 29 14:39:35 2022 .. Created module 49 | logLik: -80555550.5823023
Sat Jan 29 14:39:45 2022 .. Created module 50 | logLik: -80547900.5944912
Sat Jan 29 14:39:55 2022 .. Created module 51 | logLik: -80543667.0861826
Sat Jan 29 14:40:05 2022 .. Created module 52 | logLik: -80538240.9487562
Sat Jan 29 14:40:13 2022 .. Created module 53 | logLik: -80526553.0783506
Sat Jan 29 14:40:22 2022 .. Created module 54 | logLik: -80524128.2320462
Sat Jan 29 14:40:29 2022 .. Created module 55 | logLik: -80521149.8187733
Sat Jan 29 14:40:36 2022 .. Created module 56 | logLik: -80520196.8670347
Sat Jan 29 14:40:44 2022 .. Created module 57 | logLik: -80520670.0451455
Sat Jan 29 14:40:51 2022 .. Created module 58 | logLik: -80519862.6616851
Sat Jan 29 14:41:04 2022 .. Created module 59 | logLik: -80517061.981843
Sat Jan 29 14:41:13 2022 .. Created module 60 | logLik: -80512551.9269634
Sat Jan 29 14:41:23 2022 .. Created module 61 | logLik: -80511502.6443224
Sat Jan 29 14:41:33 2022 .. Created module 62 | logLik: -80508622.7458812
Sat Jan 29 14:41:40 2022 .. Created module 63 | logLik: -80510057.0595091
Sat Jan 29 14:41:53 2022 .. Created module 64 | logLik: -80512151.632463
Sat Jan 29 14:42:01 2022 .. Created module 65 | logLik: -80512798.3010048
Sat Jan 29 14:42:10 2022 .. Created module 66 | logLik: -80512476.0690908
Sat Jan 29 14:42:20 2022 .. Created module 67 | logLik: -80510657.6107395
Sat Jan 29 14:42:29 2022 .. Created module 68 | logLik: -80506678.883933
Sat Jan 29 14:42:36 2022 .. Created module 69 | logLik: -80508798.1462525
Sat Jan 29 14:42:47 2022 .. Created module 70 | logLik: -80509796.2184587
Sat Jan 29 14:43:02 2022 .. Created module 71 | logLik: -80505364.0587875
Sat Jan 29 14:43:10 2022 .. Created module 72 | logLik: -80505570.7437202
Sat Jan 29 14:43:18 2022 .. Created module 73 | logLik: -80491110.8130347
Sat Jan 29 14:43:38 2022 .. Created module 74 | logLik: -80491257.646636
Sat Jan 29 14:43:46 2022 .. Created module 75 | logLik: -80488948.8808976
Sat Jan 29 14:43:54 2022 .. Created module 76 | logLik: -80487673.7593699
Sat Jan 29 14:44:03 2022 .. Created module 77 | logLik: -80486994.0667674
Sat Jan 29 14:44:15 2022 .. Created module 78 | logLik: -80484812.8283321
Sat Jan 29 14:44:33 2022 .. Created module 79 | logLik: -80485582.00305
Sat Jan 29 14:44:45 2022 .. Created module 80 | logLik: -80486841.6288006
Sat Jan 29 14:44:53 2022 .. Created module 81 | logLik: -80484911.5126116
Sat Jan 29 14:45:02 2022 .. Created module 82 | logLik: -80483832.3752569
Sat Jan 29 14:45:11 2022 .. Created module 83 | logLik: -80486853.7127262
Sat Jan 29 14:45:23 2022 .. Created module 84 | logLik: -80488758.5844316
Sat Jan 29 14:45:31 2022 .. Created module 85 | logLik: -80489269.5300034
Sat Jan 29 14:45:42 2022 .. Created module 86 | logLik: -80490651.5550723
Sat Jan 29 14:45:53 2022 .. Created module 87 | logLik: -80492134.365323
Sat Jan 29 14:46:02 2022 .. Created module 88 | logLik: -80493649.0290316
Sat Jan 29 14:46:12 2022 .. Created module 89 | logLik: -80494326.8777699
Sat Jan 29 14:46:21 2022 .. Created module 90 | logLik: -80495058.3381338
Sat Jan 29 14:46:33 2022 .. Created module 91 | logLik: -80495247.3955497
Sat Jan 29 14:46:42 2022 .. Created module 92 | logLik: -80492982.8600317
Sat Jan 29 14:46:52 2022 .. Created module 93 | logLik: -80491967.5491405
Sat Jan 29 14:47:01 2022 .. Created module 94 | logLik: -80488915.1801846
Sat Jan 29 14:47:16 2022 .. Created module 95 | logLik: -80486765.4876996
Sat Jan 29 14:47:25 2022 .. Created module 96 | logLik: -80488505.476246
Sat Jan 29 14:47:34 2022 .. Created module 97 | logLik: -80487308.9445487
Sat Jan 29 14:47:43 2022 .. Created module 98 | logLik: -80489449.025699
Sat Jan 29 14:48:01 2022 .. Created module 99 | logLik: -80492148.7167353
Sat Jan 29 14:48:14 2022 .. Created module 100 | logLik: -80493412.0487448
Sat Jan 29 14:48:27 2022 .. Created module 101 | logLik: -80493122.8867567
Sat Jan 29 14:48:35 2022 .. Created module 102 | logLik: -80494168.1795715
Sat Jan 29 14:48:58 2022 .. Created module 103 | logLik: -80491790.2774153
Sat Jan 29 14:49:10 2022 .. Created module 104 | logLik: -80492427.987431
Sat Jan 29 14:49:21 2022 .. Created module 105 | logLik: -80491560.1555778
Sat Jan 29 14:49:35 2022 .. Created module 106 | logLik: -80490717.8243352
Sat Jan 29 14:49:46 2022 .. Created module 107 | logLik: -80492270.3434763
Sat Jan 29 14:49:55 2022 .. Created module 108 | logLik: -80495473.3107824
Sat Jan 29 14:50:17 2022 .. Created module 109 | logLik: -80496863.8633222
Sat Jan 29 14:50:31 2022 .. Created module 110 | logLik: -80494634.0799794
Sat Jan 29 14:50:42 2022 .. Created module 111 | logLik: -80497873.3084055
Sat Jan 29 14:50:52 2022 .. Created module 112 | logLik: -80497778.32996

plotGridSearchPerplexityDiff(moduleSplit, altExpName ="featureSubset", sep = 5)

L <- 125

temp <- subsetCeldaList(moduleSplit, list(L = L))
sce <- recursiveSplitCell(sce, sampleLabel = sce$mouse_id, initialK = 3, maxK = 40, yInit = celdaModules(temp))

==================================================
Starting recursive cell population splitting.
==================================================
Sat Jan 29 18:51:35 2022 .. Collapsing to 125 modules
Sat Jan 29 18:51:35 2022 .. Initializing with 3 populations
Sat Jan 29 19:04:28 2022 .. Current cell population 4 | logLik: -81853506.7819414
Sat Jan 29 19:06:43 2022 .. Current cell population 5 | logLik: -81585864.2994864
Sat Jan 29 19:07:17 2022 .. Current cell population 6 | logLik: -81334831.6896297
Sat Jan 29 19:09:55 2022 .. Current cell population 7 | logLik: -81173767.6214186
Sat Jan 29 19:11:18 2022 .. Current cell population 8 | logLik: -81035609.9914841
Sat Jan 29 19:11:50 2022 .. Current cell population 9 | logLik: -80932773.0890068
Sat Jan 29 19:12:50 2022 .. Current cell population 10 | logLik: -80825590.6792533
Sat Jan 29 19:16:39 2022 .. Current cell population 11 | logLik: -80759974.7812034
Sat Jan 29 19:17:09 2022 .. Current cell population 12 | logLik: -80705667.6693766
Sat Jan 29 19:17:40 2022 .. Current cell population 13 | logLik: -80652995.681211
Sat Jan 29 19:19:00 2022 .. Current cell population 14 | logLik: -80613063.9037478
Sat Jan 29 19:19:56 2022 .. Current cell population 15 | logLik: -80577385.7907885
Sat Jan 29 19:21:19 2022 .. Current cell population 16 | logLik: -80554725.1860722
Sat Jan 29 19:21:51 2022 .. Current cell population 17 | logLik: -80536571.2741023
Sat Jan 29 19:23:14 2022 .. Current cell population 18 | logLik: -80516481.1684355
Sat Jan 29 19:24:35 2022 .. Current cell population 19 | logLik: -80502227.5914425
Sat Jan 29 19:25:55 2022 .. Current cell population 20 | logLik: -80484716.528201
Sat Jan 29 19:26:54 2022 .. Current cell population 21 | logLik: -80471172.508403
Sat Jan 29 19:27:49 2022 .. Current cell population 22 | logLik: -80455240.7238886
Sat Jan 29 19:28:44 2022 .. Current cell population 23 | logLik: -80441789.3467488
Sat Jan 29 19:29:18 2022 .. Current cell population 24 | logLik: -80433928.1822051
Sat Jan 29 19:30:23 2022 .. Current cell population 25 | logLik: -80422869.8324057
Sat Jan 29 19:31:18 2022 .. Current cell population 26 | logLik: -80416193.3759105
Sat Jan 29 19:31:56 2022 .. Current cell population 27 | logLik: -80410207.8771751
Sat Jan 29 19:32:29 2022 .. Current cell population 28 | logLik: -80402701.7079324
Sat Jan 29 19:33:30 2022 .. Current cell population 29 | logLik: -80393957.8092313
Sat Jan 29 19:34:32 2022 .. Current cell population 30 | logLik: -80388002.9380147
Sat Jan 29 19:35:12 2022 .. Current cell population 31 | logLik: -80382757.2780166
Sat Jan 29 19:35:48 2022 .. Current cell population 32 | logLik: -80374059.6407028
Sat Jan 29 19:36:49 2022 .. Current cell population 33 | logLik: -80368504.8712624
Sat Jan 29 19:37:23 2022 .. Current cell population 34 | logLik: -80361943.8602602
Sat Jan 29 19:37:58 2022 .. Current cell population 35 | logLik: -80357549.9265348
Sat Jan 29 19:38:36 2022 .. Current cell population 36 | logLik: -80352131.0343025
Sat Jan 29 19:39:12 2022 .. Current cell population 37 | logLik: -80345876.7648695
Sat Jan 29 19:39:46 2022 .. Current cell population 38 | logLik: -80341598.7028261
Sat Jan 29 19:40:25 2022 .. Current cell population 39 | logLik: -80336731.8943073
Sat Jan 29 19:41:06 2022 .. Current cell population 40 | logLik: -80331610.8282257
Sat Jan 29 19:41:06 2022 .. Calculating perplexity
==================================================
Completed recursive cell population splitting. Total time: 53.47163 mins
==================================================

plotGridSearchPerplexityDiff(sce, altExpName = "featureSubset", sep = 5)

sce <- subsetCeldaList(sce, list(L=125, K=35))

sce <- celdaUmap(sce)
plotDimReduceCluster(sce, reducedDimName = "celda_UMAP", size = 0.25)

markers <- c("Sftpa1", "Scgb1a1", "Foxj1", "Cd79a", "Cd3e", "Marco", "S100a8", "Hba-a1",  "Col1a1", "Pecam1", "Mki67", "Nkg7", "Fcer1a", "Acta2")
plotDimReduceFeature(sce, reducedDimName = "celda_UMAP", features = markers, size = 0.25)

for(i in seq(L)) {
  print(plotDimReduceModule(sce, reducedDimName = "celda_UMAP", module = i, size = 0.25))
}

for(i in seq(L)) {
  print(moduleHeatmap(sce, featureModule = i))
}

sce <- celdaUmap(sce)
labels <- c("1: RBC",
            "2: Club",
            "3: AT2",
            "4: AT2",
            "5: Cd14/Cd16/S100a8 Monocyte", 
            "6: Myeloid", 
            "7: Myeloid", 
            "8: B-cell", 
            "9: Cd8a T-cell", 
            "10: Cd8a T-cell", 
            "11: Cytotoxic T-cell", 
            "12: NK-cell", 
            "13: Myeloid",
            "14: Dendritic",
            "15: Myeloid",
            "16: Proliferating T-cell",
            "17: Pecam1 endothelial",
            "18: Transitional cell",
            "19 Proliferating alveolar macrophage",
            "20 Alveolar macrophage",
            "21 Alveolar macrophage",
            "22 Cd14/Cd16 Myeloid",
            "23 Cd14/Cd16 Myeloid",
            "24 Cd14/Cd16 Myeloid",
            "25 Cd14/Cd16 Myeloid",
            "26 Cd14/Cd16 Myeloid",
            "27 Ciliary",
            "28 Cd14/S100a8 Monoctye",
            "29 Pecam1 endothelial",
            "30 Pecam1 endothelial",
            "31 Col1a1 fibroblast" ,
            "32 Col1a1 fibroblast",
            "33 Acta2 smooth muscle",
            "34 Col1a1 fibroblast",
            "35 Col1a1 fibroblast",)
g <- plotDimReduceCluster(sce, reducedDimName = "celda_UMAP", size = 0.25)
library(ggplot2)
g <- g + scale_color_manual(labels = labels,
    values = distinctColors(length(labels)))
g